Evolving a puncture black hole with fixed mesh refinement 
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We present an algorithm for treating mesh refinement interfaces in numerical relativity. We discuss 
the behavior of the solution near such interfaces located in the strong field regions of dynamical 
black hole spacetimes, with particular attention to the convergence properties of the simulations. 
In our applications of this technique to the evolution of puncture initial data with vanishing shift, 
we demonstrate that it is possible to simultaneously maintain second order convergence near the 
puncture and extend the outer boundary beyond lOOM, thereby approaching the asymptotically flat 
region in which boundary condition problems are less difficult and wave extraction is meaningful. 



I. INTRODUCTION 

Numerical relativity, which comprises the solution of 
Einstein's equations on a computer, is an essential tool 
for understanding the behavior of strongly nonlinear dy- 
namical gravitational fields. Current grid-based formula- 
tions of numerical relativity feature 17 or more coupled 
nonlinear partial differential equations that are solved us- 
ing finite differences in 3 spatial dimensions (3-D) plus 
time. The physical systems described by these equations 
generally have a wide range of length and time scales, 
and realistic simulations are expected to require the use 
of some type of adaptive gridding in the spacetime do- 
main. 

A primary example of the type of physical system to 
be studied using numerical relativity is the final merger 
of two inspiraling black holes, which is expected to be a 
strong source of gravitational radiation for ground-based 
detectors such as LIGO and VIRGO, as well as the space- 
based LISA 01 . The individual black hole masses Mi 
and M2 set the scales for the binary in the source interac- 
tion region, and we can expect both spatial and temporal 
changes on these scales as the system evolves. The binary 
must be evolved for a time t ~ lOOOM, M ~ Mi -|- M2, 
starting from an orbital separation ~ lOM to simulate 
its final few orbits followed by the plunge and ringdown. 
This orbital region is surrounded by the wave zone with 
features of scale ~ lOOM, where the outgoing signals take 
on a wave-like character and can be measured. Accom- 
plishing realistic simulations of binary black hole mergers 
on even the most powerful computers clearly requires the 
use of variable mesh sizes over the spatial grid. 

Adaptive mesh refinement (AMR) was first applied in 
numerical relativity to study critical phenomena in scalar 
field collapse in 1-D j2| ; several other related studies have 
also used AMR, most recently in 2-D Q. AMR has also 
been used in 2-D to study the evolution of inhomogeneous 
cosmologies 0, @- In the area of black hole evolution. 



AMR was first applied to a simulation of a Schwarzschild 
black hole 0. Fixed mesh refinement (FMR) was used 
to evolve a short part of a (nonequal mass) binary black 
hole merger , an excised Schwarzschild black hole in an 
evolving gauge 0, and orbiting, equal mass black holes 
in a co-rotating gauge . AMR has also been used to set 
binary black hole initial data [Tol[Tl| . The propagation of 
gravitational waves through spacetime has been carried 
out using AMR, first using a single 3-D model equation 
describing perturbations of a Schwarzschild black hole 
[T^ and later in the 3-D Einstein equations Grav- 
itational waves have also been propagated across fixed 
mesh refinement boundaries, with a focus on the inter- 
polation conditions needed at the mesh boundaries to 
inhibit spurious reflected waves [T^ . 

Realistic simulations of the final merger of binary black 
holes are likely to require a hierarchy of grids, using both 
FMR and AMR. The source region would have the finest 
grids, and would be surrounded by successively coarser 
grids, encompassing the orbital region and extending into 
the wave zone out to distances > lOOM. Evolving dy- 
namical gravitational fields using such a mesh refinement 
hierarchy poses a number of technical challenges. For ex- 
ample, the gravitational waves produced by the sources 
will originate as signals in the near zone and need to cross 
fixed mesh refinement boundaries to reach the wave zone. 
In addition, Coulombic-like signals that may vary with 
time but are not wavelike in character, such as are pro- 
duced by the gravitational potential around black holes, 
can stretch across mesh boundaries. Inappropriate inter- 
polation conditions at refinement boundaries can lead to 
spurious reflection of signals at these interfaces; c/. [l^ . 
Additional complications can arise when the grid refine- 
ment is adaptive. 

In this paper, we use the evolution of a single 
Schwarzschild black hole with FMR as a numerical labo- 
ratory. We represent the black hole as a puncture with- 
out excision and use gauges with zero shift in which the 
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solution undergoes significant evolution. This tests the 
ability of our code to handle dynamically changing space- 
times in the vicinity of mesh refinement boundaries. Us- 
ing a hierarchy of fixed mesh refinements, we are able 
to resolve the strong field region near the puncture (and 
demonstrate the convergence of the solution in this re- 
gion) while locating the outer boundary at > lOOA/. In 
Sec. ^ wc describe our methodology, including the nu- 
merical implementation. The treatment of mesh refine- 
ment boundaries is discussed in Sec. Illll Black hole evo- 
lutions with FMR are presented in Sec. II VI examples are 
given of evolutions using geodesic slicing, and 1 + log 
slicing with zero shift. We conclude with a summary in 
Sec. El 



II. METHODOLOGY 



and shift /3' specify the gauge, the full derivative notation 
d/dt = d/dt — is a partial with respect to time minus 
a Lie derivative, and the notation "TF" indicates the 
trace- free part of the expression in parentheses. These 
quantities are analytically subject to the conditions 
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known respectively as the Hamiltonian and momentum 
constraints. When evaluating Eqs. (PJ we recompute 
the physical quantities from the evolved quantities using 
Eqs. O. 

Note that the covariant derivatives of the lapse in 
Eqs. Ij2bp and (|2d|) are with respect to the physical met- 
ric, and are used here for compactness. In the code, this 
is computed according to 



A. Basic Equations 

We use the BSSN form of the ADM equations [Illll. 
These equations evolve the quantities 
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written here in terms of the physical, spatial 3-metric 
7y and extrinsic curvature Kij |l7j . where all indices 
range from 1 to 3. In Eq. flejl . F^^ is the Christoffel 
symbol associated with the conformal metric 7^ j . These 
quantities evolve according to 
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+ 7'X.fc + ^f'/?'fc,-2i''^«,a (2e) 



where here and henceforth the indices of conformal quan- 
tities are raised with the conformal metric. The lapse a 
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using only the conformal BSSN quantities, and the index 
of the covariant derivative is raised on the right hand side 
of Eq. (|2b|) with the physical metric. The Ricci tensor in 
Eq. (j2d|) is also with respect to the physical metric. We 
compute it according to the decomposition 
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giving the conformal and remaining pieces of the physi- 
cal Ricci tensor. The notation Vi denotes the covariant 
derivative associated with the conformal metric. 

There are many rules of thumb in the community re- 
garding how to incorporate the constraints into the evo- 
lution equations, and, in particular, when to use the in- 
dependently evolved F* as opposed to recomputing the 
equivalent quantity from the evolved metric. We have 
made our choices manifest in the writing of the equations 
here; we largely follow the rules set out in [l^ . 



B. Numerical Implementation 

For the spatial discretization of Eqs. Q, we take the 
data to be defined at the centers of the spatial grid cells 
and use standard 0(Aa;)^ centered spatial differences 
|l9j . To advance this system in time, we use the iter- 
ated Crank-Nicholson (ICN) method with 2 iterations 
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[20j . which gives 0(At)^ accuracy. We employ interpo- 
lated Sommerfeld outgoing wave conditions at the outer 
boundary on all variables, except for the f * which 
are kept fixed at the outer boundary. Overall, the code 
is second-order convergent; specific examples of this are 
given in Sec. If VI below. 

We explicitly enforce the algebraic constraints that Aij 
is trace-free and that 7 = f after each ICN iteration. We 
enforce the trace-free condition by replacing the evolved 
variable with 

and we enforce the unit determinant condition by replac- 
ing the evolved metric with 

ii] ^ I 113 ■ 

Both of these constraints are enforced in all of the runs 
presented below. Since the F* are evolved as independent 
quantities, Eq. (|le|) acts as a further constraint on this 
system of equations. We monitor the behavior of this 
so-called constraint along with the Hamiltonian and 
momentum constraints, Eqs. (|3a|l and jSEJ 

We use the Paramesh package |^ l22f to implement 
both mesh refinement and parallelization in our code. 
Paramesh works on logically Cartesian, or structured, 
grids and carries out the mesh refinement on grid blocks. 
When refinement is needed, the grid blocks needing re- 
finement are bisected in each coordinate direction, similar 
to the technique of Ref . . 

All grid blocks have the same logical structure, with 
zones in the x-direction, and similarly for riy and 
Uz- Thus, refinement of a 3-D block produces eight child 
blocks, each having HxriyUz zones but with zone sizes in 
each direction a factor of two smaller than in the par- 
ent block. Refinement can continue as needed on the 
child blocks, with the restriction that the grid spacing 
can change only by a factor of two, or one refinement 
level, at any location in the spatial domain. Each grid 
block is surrounded by a number of guard cell layers that 
are used to calculate finite difference spatial derivatives 
near the block's boundary. These guard cells are filled 
using data from the interior cells of the given block and 
the adjacent block; see Sec, liiil 

Paramesh can be used in applications requiring FMR, 
AMR, or a combination of these. The package takes care 
of creating the grid blocks, as well as building and main- 
taining the data structures needed to track the spatial 
relationships between blocks. Paramesh handles all inter- 
block communications and keeps track of physical bound- 
aries on which particular conditions are set, guarantee- 
ing that the child blocks inherit this information from 
the parent blocks. In a parallel environment, Paramesh 
distributes the blocks among the available processors to 
achieve load balance, minimize inter-processor communi- 
cations, and maximize block locality. 

This scheme provides excellent computational scala- 
bility. Equipped with Paramesh, the scalability of our 



code has been tested for up to 256 processors and has 
demonstrated a consistently good scaling factor for both 
unigrid (uniform grid) and FMR runs. For unigrid runs, 
we started with a uniform Cartesian grid of a certain 
number of grid cells, a fixed number of timesteps, and a 
certain number of PEs (Processing Elements), and then 
increased the number of PEs to run a larger job while 
the number of grid cells per PE remained constant. In 
this situation, we expect that the total time taken to run 
the code, including the CPU time used by all of the PEs, 
should scale linearly with the size of the problem under 
perfect conditions. In reality, communication overhead 
makes the scalability less than perfect. We define the 
scaling factor to be the time expected with perfect scal- 
ing divided by the actual time taken. Using FMR, we ran 
the same simulations as in the unigrid case except that 
a quarter of the computational domain was covered by a 
mesh with twice the resolution. Despite the more com- 
plicated communication patterns, scalability in the FMR 
runs is comparable to that in the unigrid runs. The scal- 
ing factor of our code is 0.92 for unigrid runs and 0.90 
for FMR runs. 

For the work described in this paper, we are using 
FMR. For simplicity, we use the same timestep, chosen 
for stability on the finest grid and with a Courant factor 
of 0.25, over the entire computational domain; cf. [23| . 
At the mesh refinement boundaries, we use a single layer 
of guard cells. Special attention is paid to the restriction 
(transfer of data from fine to coarse grids) and prolon- 
gation (coarse to fine) operations used to set the data in 
these guard cells, as discussed in the next section. 



III. TREATMENT OF REFINEMENT 
BOUNDARIES 

Careful treatment of guard cells at mesh refinement 
boundaries is needed to produce accurate and robust nu- 
merical simulations. The current version of our code uses 
a third order^ guard cell filling scheme that is now in- 
cluded with the standard Paramesh package. This guard 
cell filling proceeds in three steps. 

The first step is a restriction operation in which interior 
fine grid cells are used to fill the interior grid cells of the 
underlying "parent" grid. The parent grid is a grid that 
covers the same domain as the fine grid but has twice 
the grid spacing. The restriction operation is depicted 
for the case of two spatial dimensions in the left panel of 



In our terminology the "order of accuracy" refers to the order 
of errors in the grid spacing. Thus, third order accuracy for 
guard cell filling means that the guard cell values have errors of 
order Ax^, where Ax is the (fine) grid spacing. (Note that third 
order accurate gu ard cell filling was termed "quadratic" guard 
cell filling in Ref. Second order accuracy for the evolution 

code means that, after a finite evolution time, the field variables 
have errors of order Ax^. 
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FIG. 1: Guard cell filling in two spatial dimensions. In these pictures, the thick vertical line represents a refinement boundary 
separating fine and coarse grid regions. The picture on the left shows the first step, in which one of the parent grid cells (grey 
square) is filled using quadratic interpolation across nine interior fine grid cells (black circles). The other parent grid cells are 
filled using corresponding stencils of nine interior fine grid cells. (The asymmetry in the left panel is drawn with the assumption 
that the fine block's center is toward the top-left of the panel.) The picture on the right shows the second step in which two 
fine grid guard cells (grey circles) are filled using quadratic interpolation across nine parent grid values (squares). These parent 
grid values include one layer of guard cells (black squares) obtained from the coarse grid region to the right of the interface, 
and two layers of interior cells (grey squares). The final step in guard cell filling (not shown in this figure) is to use "derivative 
matching" to fill the guard cells for the coarse grid. 



Fig.d 

The restriction proceeds as a succession of one- 
dimensional quadratic interpolations, and is accurate to 
third order in the grid spacing. Note that the 3-cell-wide 
fine grid stencil used for this step (nine black circles in 
the figure) cannot be centered on the parent cell (grey 
square). In each dimension the stencil includes two fine 
grid cells on one side of the parent cell and one fine grid 
cell on the other. The stencil is always positioned so that 
its center is shifted toward the center of the block (as- 
sumed in the figure to be toward the upper left). This 
ensures that only interior fine grid points, and no fine 
grid guard cells, are used in this first step. 

For the second step, the fine grid guard cells arc filled 
by prolongation from the parent grid. Before the pro- 
longation, the parent grid gets its own guard cells (black 
squares in the right panel of Fig. ^ from the neighbor- 
ing grids of the same refinement level, in this case from 
the coarse grid. The stencil used in the prolongation 
operation is shown in the right panel of Fig. The 
prolongation operation proceeds as a succession of one- 
dimensional quadratic interpolations, and is third order 
accurate. In this case, the parent grid stencil includes 
a layer of guard cells (black squares), as well as its own 
interior grid points (grey squares). At the end of this 
second step the fine grid guard cells are filled to third 
order accuracy. 



ing" at the interface.^ With derivative matching the 
coarse grid guard cell values are computed so that the 
first derivatives at the interface, as computed on the 
coarse grid, match the first derivatives at the interface 
as computed on the fine grid. The first derivative on 
the coarse grid is obtained from standard second order 
differencing using a guard cell and its neighbor across 
the interface. The first derivatives on the fine grid are 
computed using guard cells and their neighbors across 
the interface, appropriately averaged to align with the 
coarse grid cell centers. This third step fills the coarse 
grid guard cells to third order accuracy. 

An alternative to derivative matching, which we do not 
use, is to fill the coarse grid guard cells from the first layer 
of interior cells of the parent grid. However, we find that 
such a scheme leads to unacceptably large refiection and 
transmission errors for waves passing through the inter- 
face. These errors are suppressed by derivative matching. 

Why should third order guard cell filling be adequate 
to maintain overall second order accuracy? This is a 
nontrivial question, and there are certain subtleties that 
arise in our black hole evolutions. In Appendix ^ we 
present a detailed error analysis for our guard cell fill- 
ing algorithm based on simplified model equations for a 
scalar field in 1-D. This toy model shares many of the 
features of the full BSSN system and provides a useful 
guide to understanding the behavior of our black hole 
evolutions. We demonstrate that, with this algorithm. 



The third step in guard cell filling is "derivative match- 



^ In Ref. this process is referred to as "flux matching." 
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second spatial derivatives of the BSSN variables defined 
by Eqs. ^ acquire first order errors at grid points ad- 
jacent to mesh refinement boundaries. These first order 
errors show up as spikes in a convergence plot for quan- 
tities that depend on second spatial derivatives, such as 
the Hamiltonian constraint. The key result of this analy- 
sis, however, is a demonstration that the first order errors 
in second derivatives do not spoil the overall second or- 
der convergence of the evolved variables in Eqs. in 
spite of the fact that second spatial derivatives appear 
on the right-hand sides of the evolution equations |(2J); c/. 



IV. BLACK HOLE EVOLUTIONS 




Black hole spacetimes are a particularly challenging 
subject for numerical study. Astrophysical applications 
will require that our FMR implementation perform ro- 
bustly under the adverse conditions which arise in black 
hole simulations, such as strong time-dependent poten- 
tials and propagating signals that become gravitational 
waves. In this section we demonstrate that our tech- 
niques perform convergently and accurately in the pres- 
ence of strong field dynamics and singular "punctures" 
associated with black hole evolutions, and that these 
methods can be stable on the timescales required for in- 
teresting simulations. 

The puncture approach to black hole spacetimes gener- 
alizes the Brill-Lindquist prescription of initial data 
for black holes at rest. In this approach, the spacetimc 
is sliced in such a way as to avoid intersecting the black 
hole singularity, and the spatial slices are topologically 
isomorphic to R'^ minus one point, a puncture, for each 
hole. The punctures represent an inner asymptotic re- 
gion of the slice which can be conformally transformed 
to data which are regular on R'^. In this way a resting 
black hole of mass M located at r = is expressed in 
isotropic spatial coordinates by jij = ^%]^6ij, with con- 
formal factor '^BL = 1 + M/2r and Kij = 0. A direct 
generalization of this expression for the conformal factor 
can be used to represent multiple black hole punctures, 
and data for spinning and moving black holes can be con- 
structed according to the Bowen-York prescription. 

A key characteristic which makes this representation 
appropriate for spacetime simulations is that with suit- 
able conditions on the regularity of the lapse and shift, 
the evolution equations imply that time derivatives of 
the data at the puncture arc regular everywhere despite 
the blow up in "^^bl at the puncture. Numerically, we 
treat the punctures by a prescription similar to that given 
in In the BSSN formulation, this amounts to a 

splitting of the conformal factor exp(4(/)) into a regular 
part, exp(40r), and non-evolving singular part, exp(4(/)s), 
given by The numerical grid is staggered to make 

sure the puncture does not fall directly on a grid point, 
and to avoid the large finite differencing error, the deriva- 
tives of 0s arc specified analytically. 



FIG. 2: Evolution of the conformal metric component 'y^x, 
for a geodesically sliced puncture, shown at t = 0.5, 1.0, 1.5, 
2.0, and 2.5M. 



We study two test problems, each representing a 
Schwarzschild black hole in a different coordinate sys- 
tem. Both problems test the performance of our FMR 
interfaces under the condition that strong-field spacetime 
features pass through the interfaces. The first case, de- 

31 



scribed in Sec. IIV Al is a black hole in geodesic coordi- 
nates, in which the data evolve as the slice quickly ad- 
vances into the singularity. In Appendix ^ we present 
an analytic solution for the development of this space- 
time with which we can compare for a direct test of the 
simulation. Our next test, in Sec. IIVBI uses a variant 
of the "1+log" slicing condition to define the lapse, a, 
with vanishing shift. /?* = 0. This gauge choice allows 
the slice to avoid running into the singularity, but causes 
the black hole to appear to grow in coordinate space so 
that the horizon passes though our FMR interfaces. 



A. Geodesic Slicing 

We begin with the numerical evolution of a single punc- 
ture black hole using geodesic slicing. As explained in 
Appendix IbI at i = ttM the slice on which our data 
resides will reach the physical singularity"^; because we 
are not performing any excision, this sets the maximum 
duration of our evolution. Nevertheless, t sa 3M is long 
enough to test the relevant features of FMR evolution 
and provide us with a simple analytical solution. 

In these simulations we locate the puncture black hole 
at the origin. The spherical symmetry of the problem 
allows us to restrict the simulation to an octant do- 



^ In Appendix!^ we use r as the time coordinate. When using the 
results derived there in the main body of this paper, we relabel 
T — > t to conform to the notation set out in Sec. [H] 
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FIG. 3: Convergence of the errors (numerical values minus analytic values) in jxx, Axx, the Hamiltonian constraint H , and the 
momentum constraint , for a geodesically sliced puncture along the a;-axis, all at the time t = 2.5M. The solid line shows 
the errors for the highest resolution run. The errors for the medium resolution run (dashed line) and the lowest resolution run 
(dotted line) have been divided by factors of 4 and 16, respectively, to demonstrate second order convergence. Note that the 
full domain of the simulation extends to 128M. 



main with symmetry boundary conditions, thereby sav- 
ing memory and computational time. The outer bound- 
aries are planes of constant x, y, or z at 128M each. As 
noted before, one of the important uses of FMR is to en- 
able the outer boundary to be very far from the origin, 
at > lOOM. In this case we can apply our exact solution 
as the outer boundary condition, though any numerical 
effects produced at this distant boundary arc completely 
irrelevant for the most interesting strong-field region. 

In this test we are mainly interested in how the FMR 
boundaries behave near the puncture and under strong 
gravitational fields. Even with the outer boundary far 
away we can, by applying multiple nested refinement re- 
gions, highly resolve the region near the puncture, as 
is required to demonstrate numerical convergence. To 
achieve the desired resolution near the puncture, we use 8 
cubical refinement levels, locating the refinement bound- 
aries at the planes 64M, 32M, 16M, 8M, 4M, 2M and 
IM in the x-, y-, and z-directions. 

To test convergence we will examine the results of three 



runs with identical FMR grid structures, but different 
resolutions. The lowest resolution run has gridpoints 
Ax/ = Ay/ = Azf = M/IQ apart in the finest refine- 
ment region near the puncture. The medium resolution 
run has double the resolution of the first run in each re- 
finement region. The highest resolution run has twice the 
resolution of the medium resolution run, for a maximum 
resolution of Ax/ — A//64 near the puncture. The mem- 
ory demand and computational load per timestep for the 
low, medium, and high resolution runs are similar to un- 
igrid runs of 32"^, 64'^ and 128'^ gridpoints. 

Since the data in our simulations are defined at the cen- 
ters of the spatial grid cells (see Fig.QJ, we must interpo- 
late when extracting data on cuts through the simulation 
volume. We use cubic interpolation, which is accurate to 
order Ax*, to insure that the interpolation errors are 
smaller than the largest differencing errors of order Ax^ 
expected in the simulations; c/. the discussion on post- 
processing in |29l| . When interpolating at a location near 
a refinement boundary, we adjust the stencil so that the 
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interpolation involves only data points at the same level 
of refinement while still maintaining order Ax^ errors. 

In Fig. |2]we plot the conformal metric component ^xx 
for the highest resolution run along the x-axis at times 
t = 0.5, 1.0, 1.5, 2.0, and 2.5Af. Note that in these 
coordinates the event horizon is at r = 0.5Af ai t — Q 
and moving outward toward larger values of the radial 
coordinate. By t = ttM, the singularity is at coordinate 
position r = 0.5A/, so the mesh refinement interface at 
X = IM is truly in the strong field regime. Because the 
slice will hit the singularity at coordinate position x = 
0.5M, the metric grows sharply there as the simulation 
time advances. 

In the present context though, we are not so much in- 
terested in the field values of this well-studied spacetime, 
as in the simulation errors, that is, the differences be- 
tween the analytical solution presented in Appendix 
and the numerical results. These differences allow us to 
directly measure the errors in our numerical simulation. 
At late times, these errors are, not surprisingly, domi- 
nated by finite differencing errors in the vicinity of the 
developing singularity. 

The plots in Fig. |21 compare these errors along the x- 
axis at t = l.hM for the three different resolution runs 
described above, demonstrating the convergence of "^xx-, 
Axx, the Hamiltonian constraint, and the a;-component 
of the momentum constraint. In each panel, the solid 
line shows the errors for the high resolution run. The er- 
rors for the medium (dashed line) and low (dotted line) 
resolution runs have been divided by 4 and 16, respec- 
tively. That the curves shown lie nearly atop one an- 
other is an indication of second order convergence, i.e. 
that the lowest order error term depends quadratically 
on the gridspacing. Ax. That the remaining difference 
between the adjusted curves near x = Q.5AI seems also 
to decrease quadratically is an indication that the next 
significant error term is of order Ax^. We achieve conver- 
gence to the analytic solution everywhere, from very near 
X = 0, through the peak region which is approaching the 
singularity, and in the weak field region. Animations of 
these results can be found in the APS auxiliary archive 
EPAPS: see Fig. 3A and the associated animation file in 
Ref. [33. 

Our particular interest is in the region near the re- 
finement interfaces. Fig. 0] shows a close-up of ^xx near 
the refinement boundary at x = 2M. In this figure wc 
have included the values of the guardccUs used for defin- 
ing finite-difference stencils near the interface. Again, 
the curves lie nearly atop one another, indicating second 
order convergence. 

A similar close-up of the Hamiltonian constraint H re- 
quires a little more explanation. Eq. H3a|) involves up to 
second derivatives of the BSSN data, implemented by fi- 
nite differencing. Fig.|31shows the computed values of the 
Hamiltonian constraint in the vicinity of the refinement 
boundary at a; = 2M. While the result is second-order 
convergent at any specific physical point in the neigh- 
borhood of the boundary, the figure indicates that the 
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FIG. 4: Close-up of the convergence of the error (numerical 
value minus analytic value) in "yxx J for a geodesically sliced 
puncture, &t t — 2.5M, in the vicinity of the refinement 
boundary at a; = 2M. The errors for the high resolution 
run are shown by the solid line; the errors for the medium 
(dashed line) and low (dotted line) resolution runs have been 
divided by factors of 4 and 16, respectively. In this plot, we 
also show the location of the data points, including guardcells, 
using filled circles. 
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FIG. 5: A view of the Hamiltonian constraint in the neigh- 
borhood of the FMR interface at a:: = 2M. The errors for 
the high resolution run are shown by the solid line; the errors 
for the medium (dashed line) and low (dotted line) resolution 
runs have been divided by factors of 4 and 16, respectively. 
As discussed in the text, data points nearest to the interface 
converge at one order lower than in the rest of the domain. 



sequence of values computed at the nearest point ap- 
proaching the interface as Ax — > approaches zero at 
only first order. This is as expected according to the dis- 
cussion in ApDcndixfXl fcf. Fig. ll2|l for a derived quantity 
involving second derivatives. We have specifically verified 
that, as with jxx, all BSSN variables converge to second 
order at the refinement boundaries. 

We also examined the simulation data along cuts away 
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from the cc-axis and have found them to be quahtativcly 
similar to those on the axis. In particular, plots and an- 
imations of the errors along the line y = z = 0.25M can 
be found in the EPAPS supplement; see Fig. 3B and the 
associated animation file in Ref. [sOI ■ This particular 1-D 
cut is instructive since it includes the strong field region 
yet has no particular symmetric relation to the solution. 
The fact that the errors along this line are qualitatively 
similar to those along the x-axis gives us confidence that 
the results we display in Fig. El are not subject to acci- 
dental cancellations due to octant symmetry boundary 
conditions that might produce artificially small errors. 

We have also examined the Li and L2 norms of the er- 
rors in basic variables and constraints to assess the overall 
properties of the simulation. Representative results are 
shown in Fig. El where the top panel displays the conver- 
gence behavior of the L2 norm of of the error in "f^x and 
the bottom panel the convergence of the L2 norm of H. 
The errors for the medium (dashed line) and low (dot- 
ted line) resolution runs have been divided by 4 and 16, 
respectively. These curves lie nearly atop the errors for 
the high resolution run (solid line), indicating the second 
order convergence of these error norms; see Appendix (21 
and Eqs. ^ and (|^ . 

B. 1 + log slicing 

Having rigorously tested the code against an analytic 
solution, we now use a different coordinate condition to 
study a longer-lived run with nontrivial, nonlinear dy- 
namical behavior in the region of FMR interface bound- 
aries. For this purpose, we again use zero shift but with 
a modified "l-f-log" slicing condition given by 

^ = -2«*|^X, (8) 

where insertion of the factor ^'bl = l + M/2r, originally 
recommended by [l^ , has proven to enhance convergence 
near the puncture in our simulations. For the numerical 
experiments with 1-flog slicing, the grid structure, in- 
cluding the locations of the mesh refinement interfaces, 
is the same as in the geodesic slicing case. We carry out 
three runs, with low, medium and high resolution defined 
as before. 

A l-l-log evolution serves as an excellent numerical ex- 
periment to test the robustness of our mesh refinement 
interfaces. The 1-t-log family has been well-studied in un- 
igrid runs in the past, so the generic behavior of this co- 
ordinate system is known and provides a general context 
for comparison with our mesh refinement results. Be- 
cause the l-|-log slicing is singularity avoiding, in contrast 
to the geodesic slicing case, simulations in a H-log gauge 
are known to last ~ 30Af — 40M, giving us an opportunity 
to study the properties of our mesh refinement interfaces 
in longer duration runs. Finally, as shown by Fig. 
as the lapse (right panel of the figure) collapses around 
the singularity, a strong gradient region in the metric 




Time (M) 

FIG. 6: Convergence of the L2 norms of the errors in ^x:^ 
and the Hamiltonian constraint H for a geodesically sliced 
puncture. The numerical data points are here marked with 
the symbol "x" . The solid lines show the errors for the highest 
resolution run. The errors for the medium resolution run 
(dashed lines) and the lowest resolution run (dotted lines) 
have been divided by factors of 4 and 16, respectively, to 
demonstrate second order convergence. 



(left panel of the figure) moves outward, passing through 
mesh refinement boundaries in the process. According to 
unigrid runs already in the literature {e.g., |l8l |'l choosing 
an appropriate shift, such as the Gamma-driver shift, 
would cause the evolution to freeze, preventing catas- 
trophic growth in the metric functions and confining the 
strong field behavior to the region r < lOM. This also 
increases the stable evolution time of the simulations. 
For our purposes here, however, we choose to let the 
strong gradient region move outward because we specif- 
ically wish to study how well the mesh refinement inter- 
faces handle a strong dynamical potential on timescales 
t > lOM. We consider this an important test, since such 
phenomena may develop near refinement boundaries in 
the course of realistic astrophysical simulations of multi- 
ple black holes. 

Having made this choice, we expect to see exactly what 
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FIG. 7: A time evolution sequence for the conformal metric 'y^x and the lapse a using the variant of 1+log slicing given in 
Eq. 0. Results are shown for the highest resolution run. 



appears in Fig. |7| The metric function '^xx (left panel) 
grows due to well-understood grid stretching related to 
the collapse of the lapse (right panel) and the fact that 
grid points arc falling into the black hole. The peak of 
the metric simultaneously moves to larger coordinate po- 
sition. We expect, therefore, that at some point certain 
regions of the simulations will no longer exhibit second 
order convergence because the gradients in the metric 
simply grow too large, because the peak of the metric 
moves into a region of lower refinement that cannot re- 
solve the gradients already present in the metric at that 
point, or because of a combination of the two. The sim- 
ulations in this gauge, nonetheless, remain second order 
convergent long enough for us to study the effects of the 
strong potential passing through the innermost mesh re- 
finement interfaces. 

Because we do not have an analytic solution for the 
l-|-log case to use in our convergence tests, we show three- 
point convergence plots instead. Specifically, for a given 
field /, we plot [f\ow ~ /mod)/4 using a dashed Hnc and 
/med — /high usiug a solid liuc. Since the three different 
resolutions "low" , "medium" , and "high" are related to 
each other by factors of two, the two lines in each panel 
should overlay exactly for perfect second order conver- 
gence. 

Fig.|Slshows such a three-point convergence plot for 
and Axx for a 1-D cut along the x-axis. The left panels, 
showing data from t = 8M, demonstrate that the metric 
and other variables are second order convergent every- 
where at that time. Overall, we continue to see second 
order convergence in the evolved variables, constraints, 
and norms until t ^ lOM. 

The convergent behavior starts to break down around 
t ~ lOAf due to difficulties with resolving the sharp fea- 
ture in the metric. In the region IM < a; < 2M, be- 
tween the first and second FMR boundaries, the peak 
itself grows sharply and the coarser grid is not sufficient 
to provide the resolution needed for convergent behavior. 



For 2M < a; < 4Af , the grid is again coarsened by a fac- 
tor of 2 and is not able to resolve adequately the steep 
gradient on the leading edge of the metric peak. A snap- 
shot at t — 16M is shown in the right panels of Fig. |H1 
by this time, the peak of the metric has passed through 
two refinement interfaces (at x = \M and x = 2M). The 
time development of these errors, and in particular their 
departure from second order convergence, can be seen in 
the animations available in the EPAPS supplement: see 
Fig. 8 A and the associated animation file in Ref . [Sfl . 

Throughout the duration of the runs the region .t > 5 
does remain second order convergent, even though the 
grid is further coarsened by factors of two at x = 8M, 
16M, 32M, and 64M, since all the fields change very 
slowly as they approach the asymptotically flat regime. 
The simulations will continue to run stably past this 
point (to approximately t w 35M), but the resolution in 
the regions to the right of the interface at x = \M is not 
sufficient to produce convergent results, as was expected. 

The Hamiltonian and momentum constraints along the 
X-axis are shown in Fig. |51 Three curves are plotted in 
each panel. The errors for the highest resolution run 
are given by the solid line. The errors for the medium 
(dashed line) and low (dotted line) resolution runs have 
been divided by factors of 4 and 16, respectively. The 
constraints are second order convergent in the bulk for 
times t ^ lOM, when the resolution is sufficient to han- 
dle the growing feature in the metric (left panels). As ex- 
pected, H exhibits first order convergent spikes at mesh 
refinement interfaces; c/. Fig. |S1 and Appendix ^ For 
t > lOAf , as the peak of the fields propagates into the 
coarser grid regions past x = 2M, the lowest resolutions 
are not sufficient to resolve the rising slope of the metric, 
and, like the evolved variables (Fig.lSJ, the constraints no 
longer demonstrate second order convergence. The right 
panels of Fig.|51show the constraints att = 16M, right af- 
ter the peak of the metric passes through the refinement 
interface at x = 2M. See Fig. 9A and the associated 
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FIG. 8: Three-point convergence plots for the BSSN variables 'y^x and A^x along the a;-axis in the 1+log slicing runs at times 
t = 8M (left panels) and t = 16Af (right panels). For a given field /, the dashed line shows (/low — /mcd)/4 and the solid line 

shows /med " /high- 



animation file in Ref. jSOj for animations of these data. 

The behavior of the simulations at locations away from 
the X-axis is qualitatively similar to that shown in Figs.|Sl 
and|5| Plots and animations of the errors along the line 
y = z = Q.25M are available in the EPAPS supplement; 
see Figs. 8B and 9B and their associated animation files 
in Ref. [s^. 

We have also examined the Li and L2 norms of the 
errors to assess the overall behavior of these runs, and 
display representative results in Fig. 1101 The L2 norms of 
the errors in the basic variables ^xx and A^x are shown in 
the left top and bottom panels, respectively, using 3-point 
convergence plots. The dashed lines show the difference 
between the low and medium resolution results divided 
by 4, and the solid lines show the difference between the 
medium and high resolution results, demonstrating the 
overall second order convergence of these simulations at 
early times. The Li norm of H is displayed in the top 
right panel, where the solid line gives the errors for the 
high resolution run. The errors for the medium (dashed 
line) and low (dotted line) resolution runs have been di- 
vided by factors of 4 and 16, respectively to show second 



order convergence, as expected from Eq. ljC7p . In the 
lower right panel the L2 norm of H is shown, with the 
solid line giving the results for the high resolution run. 
As discussed in Appendix |0 the errors for the medium 
(dashed line) and low (dotted line) resolution runs have 
been divided by factors of 2^^^ and 4'^/^ = 8 to account 
for the effects of significant first order convergent errors in 
H at the mesh refinement boundaries, in addition to the 
second order convergent errors in the bulk; see Eq. IjCSp . 

One final feature of these simulations, the high fre- 
quency noise near the origin seen in the right panels of 
Fig. |51 requires some explanation. First of all, it is not 
related to the presence of the refinement boundaries; in 
particular, we have reproduced it in unigrid runs and 
with an independently-written, 1-D (spherically symmet- 
ric) code. Higher resolution exacerbates this problem: 
both the frequency and amplitude of the noise increase 
with resolution. We have found the location of this noise 
to be independent of resolution and the number and po- 
sitions of FMR boundaries. 

This feature, which we call the "point-two A/ prob- 
lem," originates around r 0.2M. It becomes most evi- 
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FIG. 9: Convergence plots for the Hamiltonian constraint H and the momentum constraint along the a;-axis in the 1+log 
slicing runs at times t = 8M and t = 16M. The solid lines show the errors for the high resolution run. The errors for the 
medium (dashed line) and low (dotted line) resolution runs have been divided by factors of 4 and 16, respectively. 



dent at times t > lOM, first appearing in the lapse and 
K, which are directly coupled, and then eventually mix- 
ing into all of the extrinsic curvature variables. For the 
duration of the evolutions the noise remains within the 
region O.OM ^ x ^ 0.5M. Outside this region, all basic 
variables demonstrate satisfactory second order conver- 
gence, including at refinement boundaries, up to times 
t ~ lOM. 

Having chosen a generally accepted gauge, and having 
focused on effects of the mesh refinement interfaces in 
this work, we have not fully investigated the cause of 
nor possible remedies for this apparent pathology. We 
note it here, however, as an interesting topic for future 
investigation. 



V. SUMMARY 

This paper demonstrates that fixed mesh refinement 
boundaries can be located in the strong field region of 
a dynamical black hole spacetime when the interface 
conditions arc handled properly. This result was ver- 



ified through simulation of a Schwarzschild black hole 
in geodesic coordinates, for which we have an analytic 
solution for comparision, and through simulations of 
Schwarzschild in a variation of the 1-l-log (singularity 
avoiding) slicing with zero shift. Mesh refinement tech- 
nology, therefore, is a viable to way to use computational 
resources more efficiently, and to simulate the very large 
spatial domains needed to compute the dynamics of the 
source interactions and allow extraction of the resulting 
gravitational waveforms. 

Our method for handling the interface conditions, 
based in part on the Paramesh infrastructure, is detailed. 
For these simulations we find that, in handling the inter- 
face condition between FMR levels, third order guard cell 
filling is sufficient for overall second order accuracy in the 
simulations. By nesting several levels of mesh refinement 
regions, we are able to resolve the puncture convergently 
while simultaneously pushing the outer boundary of our 
domain to 128M and keeping the computational problem 
size modest. We estimate that for only a 12% increase 
in the computational size of the problem, we could push 
the outer boundary to 256Af ; moving the outer boundary 
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FIG. 10: Convergence behavior of the L\ and L2 norms of the errors for the runs with 1 + log shcing. The left panel shows 
the L2 norms of the errors in '^^x (top) and A^x (bottom), where the dashed lines show the the difference between the low 
and medium resolution results divided by 4, and the solid lines show the difference between the medium and high resolution 
results. The top right panel shows the L\ norm of H, with the solid line giving the errors for the high resolution run; the errors 
for the medium (dashed line) and low (dotted line) resolution runs have been divided by factors of 4 and 16, respectively. The 
lower right panel shows the L2 norm of H , with the solid line giving the results for the high resolution run. In this panel, the 
errors for the medium (dashed line) and low (dotted line) resolution runs have been divided by factors of 2"^''^ and 4"^''^ = 8, 
respectively. 



out even farther will be possible for production runs on 
larger machines. Combined with our earlier results show- 
ing that gravitational waves pass through such FMR in- 
terfaces without significant reflections [lj|, we have now 
studied, in detail, the effects of FMR interfaces on the 
two primary features, waves and time- varying strong po- 
tentials, of astrophysically interesting spacetimes. 

In this paper, we have evolved single black holes using 
gauges with zero shift in order to produce test problems 
in which strong-field spacctime features with steep gradi- 
ents pass through mesh refinement interfaces. In more re- 
alistic, astrophysical simulations of multiple black holes, 
we expect to use non-zero shift prescriptions. While 
a shift vector will allow us to control certain aspects 
of the dynamics, we still expect to find some strong, 
time- varying signals to propagate across mesh refinement 
boundaries. We are currently implementing non-zero 
shift conditions into our FMR evolutions and will report 



on this work in a separate publication. 
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APPENDIX A: ERROR ANALYSIS OF GUARD 
CELL FILLING SCHEME 

To help pave the way for understanding the behavior 
of our black hole evolutions near mesh refinement bound- 
aries, we provide here a detailed analysis of a toy model 
for a scalar field in one spatial dimension, using the same 
third order guard cell filling algorithm detailed in Sec. lIIII 
The model equations are 

Tp = TT (Ala) 
7T = V" , (Alb) 

where the dot denotes a time derivative and primes de- 
note space derivatives. These equations can be solved nu- 
merically using the same twice-iterated Crank-Nicholson 
algorithm used to evolve our black hole spacetimes. The 
fields at timestcp n + 1 are given in terms of the fields at 
timestep n by 

rp^+^ = + AtTT^I + —D^ij^l + —D\:^ (A2a) 

Af'^ Af^ 
ttJ+i = + AtD^i^:! + —D'^TT^ + —D*'i(;fA2h) 

where is a finite difference operator approximating 
the second spatial derivative, and Z?^ = {D'^Y . 

Consider for the moment a uniform spatial grid. If 
is the usual second order accurate centered difference 
operator, the dominant source of error for ^^^^ comes 
from the term proportional to At^. This term has the 
wrong numerical coefficient as compared to the Taylor 
series expansion of the exact solution. The dominant 
sources of error for Tr""*"^ come from the term proportional 
to A^"^, which also has the wrong numerical coefficient, 
and from the second order error in D'^ip^ . For a uniform 
grid the dominant error in D^ipj' is D'^ipJ Ax'^/12, so the 
leading errors for a single timestep are 

(V-r')-- = (A3a) 
{n^+\rr = ^AtiAf + Ax^)D*'ipJ. (A3b) 

For At ~ Ax, each of these one-time-step errors is pro- 
portional to Aa;^. If we evolve the initial data to a finite 
time T, the 0{Ax^) errors accumulate over N ^ T/At 
timcsteps resulting in second order errors."* Thus, the 



* This is a simplification. The dominant error after A'^ timesteps 
includes other terms of order Ax^ in addition to the product of 
N and the one— time— step error. These other terms include, for 
example, the product of an order A'^'^ coefficient and a one-time- 
step error of order Ax^. 



basic variables il' and n are second order convergent on a 
uniform grid. 

On a non-uniform grid, guard cell filling introduces 
errors of order Aa;'^ in ip at grid points adjacent to the 
boundary. This leads to errors of order Aa; in ZJ^i/i" 
and 1/Ax in Z3'*'0". From Eq. (jA2b|) we see that in one 
timestep tt can acquire errors of order Aa;^ . The concern 
is that these errors might accumulate over TV = T/At 
timesteps to yield first order errors. This, in fact, does 
not happen. Simple numerical experiments show that ip 
and TT are second order convergent on a non-uniform grid 
with third order guard cell filling. 

We can understand this result with the following 
heuristic reasoning. The numerical algorithm of Eq. (|A2|) 
approximates, as does any mathematically sound numeri- 
cal scheme, the exact solution of the scalar field equations 
(Eq. (|A1|) ') in which the field tt propagates along the light 
cone. The "bulk" errors displayed in Eq. l|A3b|) accumu- 
late along the past light cone to produce an overall error 
of order N Ax^ ~ Aa:^ at each spacetime point. Errors 
in guard cell filling, which occur at a fixed spatial loca- 
tion, do not accumulate over multiple timesteps since the 
past light cone of a given spacetime point will cross the 
interface (typically) no more than once. 

The characteristic fields for the system (|Alll are tt ± ■(/;' 
so that tp', like tt, propagates along the light cone. As 
a result, the value of i/j at a given spacetime point is 
determined by data from the interior of the past light 
cone. From Eq. ljA2ap we see that the one-time-step 
errors for -0 due to guard cell filling arc order Aa;'^. These 
errors can accumulate over N timcsteps to yield errors of 
order A'^ Aa;^ ^ Ax^ . 

The derivatives ?/;' and ip" are computed at finite time 
T by evolving the ip, tt system for T/At timesteps then 
taking the centered, second order accurate numerical 
derivatives of ?/'• Numerical experiments show that ip' 
and tp" , defined in this way, are second order convergent 
on a non-uniform grid with third order guard cell filling.^ 
Continuing with our heuristic discussion, we can under- 
stand the second order convergence of tp" as follows. Let 
{tpj')err ~ Ax^ dcnotc the error in ip at grid point n, 
j, where the coefficient i?" is independent of Ax. Some 
of this error is due to guard cell filling at the mesh re- 
finement interface and some is due to the accumulation 
of "bulk" errors (|A3|) . Now, the second derivative of ip, 
computed as D'^iP'} = {iP^^^ -2iP^+ip:^_^)/Ax^, wiU con- 
tain errors of the form {D'^ipJ)err = E^+i - 2£^" + 
Since the bulk errors are smooth, the bulk contribution 
to E^^^ - 2E]' + E-jLi wiU scale as Aa;^. It is also the 
case that the errors due to guard cell filling are smooth. 
This is because the value of ip at any given point is deter- 
mined by the interior of the past light cone, so its error 
includes an accumulation of guard cell filling errors along 



^ The error in ip" is fairly noisy but the overall envelope containing 
this noise is second order convergent. 
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the history of the mesh refinement interface. In the hmit 
of high resolution this accumulation of error approaches 
the same value at neighboring grid points j — 1, j, and 
j + 1. In other words, the guard cell filling contribution 
to — 2i?" + approaches zero as Ax 0. Ev- 
idently, the guard cell filling contribution, like the bulk 
contribution, scales like Ax^. 

The discussion above indicates that we can evolve the 
scalar field system of Eq. (|Aip for a finite time on a grid 
with mesh refinement, numerically compute the second 
derivative of V', and find that ■0" is second order accurate. 
Without shift, the BSSN equations IScj) and lf2(!|l arc sim- 
ilar to the scalar field equations with 7^ playing the role 
of ip and —Aij playing the role of tt. This feature was 
one of the original motivations behind the BSSN system. 
Note that the term analogous to tp" in the tt equation 
is the term ^''"^jij.im contained in the trace-free part of 
the Ricci tensor, which appears on the right-hand side 
of Eq. Ij2d|l . Obviously there are many other terms that 
appear on the right-hand side of the dAij/dt equation. 
We can model the effect of these terms by including a 
fixed function on the right-hand side of the tt equation: 

-(/) = TT (A4a) 
TT = V/' - x" • (A4b) 

We have written the fixed function as the second deriva- 
tive of X- For simplicity we choose x to depend on x 
only, the most relevant dependence for our consideration 
of behavior across spatial resolution interfaces. The gen- 
eral solution of this system is then 

^{t,x) = ij{t,x)+x{x) (A5a) 
■K{t,x) = 7f(t,a;) (A5b) 

where tt is a solution of the homogeneous wave equa- 
tion (Eq. (im ). 

The extended model system of Eq. (|A4(I can be solved 
numerically with the discretization 

= ^J] + Mn'^ + ... (A6a) 

7r;.'+i = 7r;.' + AtpV';-^'xj) + --- (A6b) 

The higher order terms in Ai, not shown here, come from 
the iterations in our iterated Crank-Nicholson algorithm. 
It is important to recognize that the x" term is expressed 
as the numerical second derivative of Xj aiid not as the 
discretization of the analytical second derivative, {x")j- 
The reason for this choice is that D^Xj mimics the effect 
of the extra terms on the right-hand side of Eq. (|2dp 
which, in our BSSN code, depend on the discrete first 
and second derivatives of the BSSN variables and F*. 

From the discussion of the wave equation (Eq. (lAip ) we 
can anticipate the results of numerical experiments with 
the model system (Eq. (|A4|) ') on a non-uniform grid. For 
arbitrary initial data -0^, tt^, the numerical solution is 
given by 

^1 = + (A7a) 
tt" = 7f" (A7b) 



where 0", tt" is the numerical solution of the homoge- 
neous wave equation (Eq. (jAip l with initial data "0° — Xj, 
TTj . The order of convergence for ip" is determined by how 
rapidly, as Ax — > 0, the numerical solution in Eq. (|A7a|) 
approaches the exact solution ip{t,x) = tp{t,x) + x(x). 
Since Xj is simply the projection of the analytic function 
x(x) onto the numerical grid, the term Xj in the numer- 
ical solution (Eq. (jATap ) does not contribute any error. 
We have already determined that on a non-uniform grid 
'0" approaches ^{t^ x) with second order accuracy. Thus, 
we expect ip^ to be second order convergent. 

What about derivatives of ■0? The order of conver- 
gence for D^ip" is found by comparing the discrete deriva- 
tive D-^ip^j = D^ipJ + D^Xj to the analytic solution 
■0' = ^' -I- Again, as we have discussed, D^ip^ ap- 
proaches ip' with second order errors. It is also easy to see 
that the numerical derivative D^Xj approaches x' with 
second order accuracy. Away from grid interfaces this is 
obviously true, assuming that is the standard second 
order accurate centered difference operator. For points 
adjacent to a grid interface, guard cell values for ip are 
filled with third order errors. These errors lead to second 
order errors in D^ip". Overall then, we expect second 
order convergence for D^ip". 

The expected convergence rates for ip and 0' are con- 
firmed by the results shown in Fig. 1111 For these numer- 
ical tests, we chose xi^) ~ exp((.T — 50)/10) and initial 
data 

0(0, x) = ioOe-(-+io)'/4oo -f e(-50)/io (A8a) 
^(0,x) = i(x + 10)e-("+i°)'/4°° . (A8b) 

Each set of curves shows the errors at three different reso- 
lutions. Ax — 5/16, 5/32, and 5/64, where Ax is the fine 
grid spacing. The evolution time is 20.83, corresponding 
to 200, 400, or 800 timesteps (depending on the resolu- 
tion) and a Courant factor of 1/3. 

The order of convergence for the second derivative of 
-0 is determined from a comparison of D'^ipJ ~ D'^ip^ + 
D^Xj and the analytic solution 0" = + x" ■ We have 
seen that D^'ipj- approaches -ip" with second order ac- 
curacy. The situation for D^Xjj however, is somewhat 
different. Away from any grid interface D'^Xj will ap- 
proach x" with second order accuracy, assuming is 
the standard second order accurate finite difference op- 
erator. But for points adjacent to the interface, and only 
those points, guard cell filling errors of order Ax'^ in x 
will lead to first order errors in D^Xj- Thus, we expect 
to find second order convergence for D'^tp'^ at all points 
except those points adjacent to the interface. Points ad- 
jacent to the interface should be first order convergent. 

Fig. El shows the results of our convergence test for 
iP'". The spikes at the interface (x = 0) appear because 
the two grid points adjacent to the interface are only first 
order convergent. Elsewhere, the plot shows second order 
convergence. 
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FIG. 11: Convergence tests for ^ and ^' . The mesh interface 
is at a; = 0, with the fine grid on the left and coarse grid on 
the right. The errors in ip and ^p' for the high resolution case 
are shown by the solid line. The errors are divided by factors 
of 4 and 16 for the middle (dashed line) and low (dotted line) 
resolution cases, respectively. 
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FIG. 12: Convergence test for -i/)". The solid curve is the 
error in ip" for the highest resolution run. As in Fig. 1111 the 
interface is at a; = and the errors for the middle (dashed 
line) and low (dotted line) resolution cases are divided by 
factors of 4 and 16, respectively. The grid points adjacent to 
the interface do not coincide because ip" is only first order 
convergent at these points. 



The behavior demonstrated in Fig. ^| also occurs in 
the BSSN system when we examine the convergence plot 
for the Hamiltonian constraint. In graphing the Hamil- 
tonian constraint we are comparing a combination 
of grid functions that includes second derivatives of the 
BSSN variables to the exact analytical solution for H, 
namely, zero. We therefore expect spikes to appear at 
interfaces in the convergence plot for the Hamiltonian 
constraint, and indeed they do (see, for example. Figs. El 
and©. 

We wish to emphasize that the lack of second order 



convergence for second spatial derivatives at grid points 
adjacent to the interfaces is not due to any error in our 
code, or shortcoming of the numerical algorithm. Since 
the undifferentiated variables are second order conver- 
gent everywhere, we can always assure second order con- 
vergence of their derivatives by using suitable finite dif- 
ference stencils. For example, in computing D'^ip" from 
■0" we can use a second order accurate one-sided operator 
that avoids using guard cell values altogether. With 
such a choice the spikes in Fig.^jdisappear, and D'^ip^j is 
everywhere second order convergent. In our BSSN code, 
it is most convenient to compute the Hamiltonian con- 
straint using the same centered difference operator 
that we use for the evolution equations. As a conse- 
quence, spikes appear at the grid interfaces in the con- 
vergence plots (Fig. 121 and O . 



APPENDIX B: ANALYTIC SOLUTION FOR 
GEODESICALLY SLICED SCHWARZSCHILD 

In a numerical simulation, geodesic coordinates are ob- 
tained by using unit lapse and vanishing shift. This im- 
plies that the grid points will follow geodesic trajecto- 
ries through the physical spacetime. We present here a 
physical derivation of the Schwarzchild spacetime met- 
ric in this well-known coordinate system, based on those 
geodesies; an alternate derivation is available in Refs. 

iiii3ii. 

The Schwarzschild geometry in standard coordinates 
is given by 



(Bl) 



where gj.j, = g^^ = (1 — 2M/R). 

To express this metric in geodesic coordinates, consider 
a spatial Cauchy surface Eq in a 4-manifold Ad and a 
congruence of radial geodesies crossing Sq. Let the affine 
parameter r for each geodesic be zero at Eg. Considering 
subsequent slices of constant proper time, we can set a 
global time t which we use to define a new foliation T,r 
of Ai. Each geodesic in the congruence is labeled by 
the coordinates of its initial "starting" point in Eg. The 
radial position p of the starting point in Eq can thus be 
promoted to a new radial coordinate on JVI to pair with 
the time coordinate t. 

Now we will derive the metric components in this r- 
p coordinate system. The affine parameter r induces 
the normalized vector n° = (9/9r)° tangent to the 
geodesic, implying that the lapse grr = ^1- Assuming 
the geodesies begin at rest so that n° is normal to Eq 
implies that g^-p = initially. Furthermore, the geodesic 
equation n'^Van^ = requires that gTp,T = 0. Thus g^-p 
(the shift) must remain zero. 

A straightforward transformation from Eq. HB1() for 
the remaining metric coefficient yields 



9pp = 



{dR/dpf 



[{dT/dT)gTTY 



(B2) 
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The term in the denominator is the energy defined for 
geodesies on this spacetime, and it is conserved along 
the geodesies: n^Vairib^'') = 0, where = {d/dTf is 
the timehke KiUing field. On Sg one can evaluate this 
energy as y^-g^r, where = gab\T=o- This gives 



dR 



dp 



2M 



dp 



A similar application of conservation of energy in 
n°-na = — 1 yields fi, 7\: 



(2M)i/2 



R / R\ R 

— 1 + arccos \ / — 

P \ Pj \ P 



= 0. (B4) 



This expression provides an implicit definition for R = 
R(p,t), which is easily inverted numerically to high pre- 
cision. 

To perform numerical evolutions the geodesic coordi- 
nates have a drawback: the physical singularity is already 
present on the initial slice (t = 0) at p = 0. We can avoid 
this problem by going to isotropic coordinates (r, 9, (j)) by 
means of the transformation p — r {1 + M /2r) . We see 
that p — > cx) both as r — > and as r oo. For real r the 
minimum value of p is p = 2M (the horizon) at r = -M/2, 
now the surface closest to the physical singularity on the 
initial shce. Substituting p = 2M into Eq. HB4|I we see 
that geodesies originating on this surface reach the phys- 
ical singularity, = 0, at time t — ttM, defining the 
maximum temporal extent of our coordinate system. 

Returning to the metric, the transformation to 
isotropic coordinates gives us 



9pp 



dR 
dp 



-] 

2r J 



(B5) 



From Eq. IjBSp we can express the final metric as: 



dR 
dp 



M 
27 



dr^ + R^dn^. (B6) 



Expressions for the extrinsic curvature, which have not 
previously appeared in the literature, can be derived in 
a similar manner. As we know from the ADM formalism 
|l7ll32l |. the extrinsic curvature can be viewed as the rate 
of change of the spatial metric 



1 dgab 



(B7) 



2 dr 

when the lapse is unity and the shift is zero. This gives 
1 



1 



d fdR^ ^ 



2r J dr \ dp 



(B8) 



To evaluate the partial derivatives in Eqs. ljB6|) 
and ljB8|) . we note that if we have a function / ~ 
f{u, v,w) — Q defining u as an implicit function of v and 



w, we can use the chain rule and the implicit function 
theorem to show that 



du 
dv 

d^u 
dwdv 



df/dv 
df/du 

dl 
du 



(B9) 



dy 

dvdw 



d^didi fdi 

du^ dv dw V du 



df\-^ / ay df d^ df 



du 



V dvdu dw ~^ dudw dv 



(BIO) 



Taking for / the left hand side of Eq. (|B4p . and, noting 
that df /dr = 1, we conclude that: 



1+^V^ 

2r J dp 



dl 
dR 



dy d^ df 



dpdR dp 



dl 
dR 



(Bll) 



and 



R 



dl 
dR 



lU^^Kge sin" (6). (B12) 



There are no off-diagonal terms. Observe that we have 
only partial derivatives of / that can be obtained analyt- 
ically from Eq. (|B4|I . and easily evaluated numerically. 



APPENDIX C: DEFINITION OF L„ NORMS AND 
SCALING PROPERTIES 

For a function / defined on a uniform grid Ax = Ay = 
Az = h, we take the i„ norm of the function to be 



l/n 



Lnlf] 



.grid 



(CI) 



where fjkm is the value of the function at grid point 
{j,k,m); cf. |33l |. If the function is defined on a non- 
uniform grid with i refinement levels, the L„ norm be- 
comes 



Lnif] = 




1/r. 



hllfj 



fcm I 



jkrn I 



grid i 



(C2) 

where hi is the cell size on the z'^ grid. In our work, the 
function / denotes an error, either derived from compar- 
ison with an analytic solution (for example, the Hamil- 
tonian constraint for all our runs, and the basic variables 
with geodesic slicing) or from comparison with a run at a 
different resolution as part of a three-point convergence 
test (for the basic variables with 1 -I- log slicing). 

It is useful to work out the scaling behavior expected 
when error norms from runs with different resolutions 
are compared. Recall that, for the runs presented in this 
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paper, hiow = 2/iniod = 4/ihigh, and the errors in basic 
variables such as "fpq and Apq are expected to scale as 
f ^ everywhere. Let N be the characteristic number 
of zones along one dimension of a simulation volume, so 
that h ~ 1/N. We focus on the Li and L2 norms, which 
are the ones generally used to examine errors in numerical 
relativity. Then 

Li-Y. ~ ^'^'^ r^f^h^ (C3) 

grid 

and 

L,^(j2h'A ^ f (C4) 

\grid / 

SO that both the Li and L2 norms should exhibit second 
order convergence in this case. Note that these expres- 
sions are valid not only for unigrid runs but also for our 
FMR simulations, since the refinement structure of the 
grid is the same in all these runs. 

This situation regarding the Li and L2 norms of the 
Hamiltonian constraint H is somewhat more compli- 
cated. As we have shown in Sec. II VI and Appendix VK\ H 
has errors that scale as f ^ h on refinement boundaries 
and as / ~ /i^ in the bulk. For the runs with geodesic 
slicing, the errors in H in the bulk near the puncture 



dominate over those at the refinement boundaries; see 
Fig. 13 Since these errors show second order convergence 
f ^ h^, wc expect that both the Li and L2 norms will 
also scale ^ h^, a,s in Eqs. (|C3p and (|C4|) . However, in 
the case of 1 -I- log slicing, the first order convergent er- 
rors on the refinement boundaries play a larger role; see 
Fig. El To account for this, we write 



\ boundary bulk / 

The number of zones on the boundary ~ N'^ while the 
number of zones in the bulk ~ N^, for sufficiently large 
N. Then 

(C6) 

This gives 

Li - (C7) 

and, since /i ^ 1, 

L2 - h^/^ (C8) 
for the scaling of the norms of £f in 1 + log slicing. 
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